Ab Initio Calculation of Finite Temperature Charmonium Potentials 
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The interquark potential in charmonium states is calculated for the first time in both the zero and 
non-zero temperature phases from a first-principles lattice QCD calculation. Simulations with two 
dynamical quark flavours were used with temperatures T in the range QAT C < T < 1.7T C , where T c 
is the deconfining temperature. The correlators of point-split operators were analysed to gain spatial 
information about the charmonium states. A method, introduced by the HAL QCD collaboration 
and based on the Schrodinger equation, was applied to obtain the interquark potential. We find a 
clear temperature dependence, with the central potential becoming flatter (more screened) as the 
temperature increases. 
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Introduction - The quark-gluon plasma (QGP) phase 
of QCD has been studied extensively both in heavy-ion 
collision experiments at RHIC [l], and the LHC |3| as 
well as in theoretical calculations. However, a complete 
understanding of this phase is still some distance away. 
Experiments are hindered by uncertainties in the phe- 
nomenology of the QGP such as the equation of state, 
transport properties, and spectral features of hadrons. 
These quantities are required to model the QGP fireball 
in heavy-ion collisions as it expands and cools back into 
the hadronic phase in order that the events in the detec- 
tors can be properly interpretted. 

One of the quantities of interest is the interquark po- 
tential in the QGP phase. A temperature dependent 
charmonium potential underlies the widely cited J/ip 
suppression model of Rcf. 4]. More recent work on 
statistical models of charmonium production [H, @ and 
studies assuming transport models of charmonium pro- 
duction @,@] lead to alternative interpretations. An anal- 
ogous suppression has recently been found in bottomo- 
nium yields in heavy-ion collisions [il [T^] ■ 

Theoretical work on the interquark potential at high 
temperature includes early models (Tl| and perturbative 
QCD calculations [12] ■ Furthermore, there have been 
some recent non-perturbative (i.e. lattice) QCD studies 
of interquark potentials which are relevant to the work 
presented here. These fall into two categories: (i) non- 
zero temperature studies of the static quark potential 
[TJ-[i3 and (ii) zero temperature studies of the potential 
between quarks with finite masses [U)]. The work pre- 
sented here is a study of the interquark potential of char- 
monium using physical charm quark masses at finite tem- 
perature and uses two flavours of light dynamical quark. 
A particular feature of our work is that our lattices are 
anisotropic which has the significant advantage that our 
correlation functions are determined at a large number 
of temporal points, hence aiding our analysis. 

The method we use is based on the HAL QCD collabo- 
ration's calculation of the internucleon potential relevant 
for nuclear physics and utilised the Schrodinger equation 



18 1 . In this work we use their "time-dependent" method 
20] to determine the real part of the interquark char- 
monium potential. In our work we do not consider the 
width of the state and therefore have access to the real 
part of the potential only. The possible limitations of 
the underlying assumption, that a nonrelativistic poten- 
tial description is valid for these temperatures and quark 
masses, is a separate issue which will not be discussed 
here. 

Our main conclusion is that the charmonium potential 
as a function of distance is steepest for low temperatures 
T, and becomes natter at large distances as T increases. 
This work extends our earlier work in (2lj . 

Time- dependent Schrodinger Equation Approach - Fol- 
lowing HAL QCD, we determine the potential using their 
"time-dependent" method [2(| . The first step is to define 
charmonium point-split operators, 



Jt(x; r) = q(x) T U{x, x + r) q(x + r), 



(1) 



where r is the displacement |3l| between the charm and 
anti-charm quark fields q and q, x is the space-time point 
(x, r) and T is a Dirac matrix used to generate vector 
(J/VO or pseudoscalar (n c ) channels. U(x, x + r) is the 
gauge connection between x and x + r. The correlation 
functions, 



C r (r,r) = ]T( J r (x, r; r) 4(0; 0)). 



(2) 



of the point-split and local operators, can be expressed 
in the usual spectral representation, 

3 J 

where the sum is over the states j with the same quan- 
tum numbers as the operator Jr, and tpj{ r ) are the cor ~ 
responding Nambu Bethe Salpeter (NBS) wavefunctions. 
N T is the number of lattice points in the temporal direc- 
tion and is related to the temperature by T = l/(a T N T ), 
where a T is the temporal lattice spacing. 
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From now on wc consider only radially symmetric (S- 
wave) states. We differentiate Eq.© w.r.t. time and ap- 
ply the Schrodinger equation which, in Euclidean space- 
time is 



V r (r) 



(4) 



where fi is the reduced mass of the cc system, fj, = \m c ~ 
jMj/,p. Ignoring the backward moving contribution, we 
obtain 



dC r (r,T) 
8t 



j 



V r (r) ) J VJT e~ EiT 



2Ei 



1 <9 2 



C r (r,r). 



(5) 



This can be trivially solved for the potential Vr{r). 

Notice that the NBS wavefunction, ip(r), is not ex- 
plicitly required in the above derivation of V(r). How- 
ever, we note that HAL QCD's original "wavefunc- 
tion" method extracts ip(r) from a fit to the large 
time behaviour of the correlation function, C(r, r) — > 
tpo(Q)ipo(r) e~ EjT , and then uses this ipo(r) as input into 
the Schrodinger equation to obtain the potential fl8j |. 
HAL QCD's time-dependent method used here has the 
distinct advantage that the correlation functions are used 
directly, without requiring a fit to the asymptotic state. 

The S-wave potential can be expressed as 



V r (r) = V c (r) + si • s 2 V s {r), 



(6) 



where Vc is the spin- independent (or "central") poten- 
tial, Vs is the spin-dependent potential, and Si j2 are the 
spins of the quarks. We have Si • s 2 = —3/4, 1/4 for the 
pscudoscalar and vector channels respectively. 

Lattice Parameters and Correlators - We performed 
lattice calculations of QCD with two dynamical flavours 
of light quark using a Wilson-type action with anisotropy 
of f = a s /a T = 6, a s ~ 0.162fm and a" 1 ~ 7.35GcV 
[22T [23l | . The other lattice parameters are listed in Table 
HI We note that the range of temperatures is from the 
confined phase up to ~ 1.7T C where T c is the deconfin- 
ing transition. The charm quark is simulated with the 
(anisotropic) clover action and its mass is set by match- 
ing the experimental rj c mass at zero temperature. 



N s N T T(MeV) T/T c N clg 



12 


80 


90 


0.42 


250 


12 


32 


230 


1.05 


1000 


12 


28 


263 


1.20 


1000 


12 


24 


306 


1.40 


500 


12 


20 


368 


1.68 


1000 



TABLE I: Lattice parameters used, including spatial and tem- 
poral dimension, N s and N T , temperature, and number of 
configurations, N c { s . 
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FIG. 1: The results for the central potential Vc(r) obtained 
from Eq. © for (a) T = 0.42T C and (b) 1.68T C . The horizon- 
tal axis is the Euclidean time, r, appearing in Eq. ((Sj 



Results - We now apply Eq. ([5]) to obtain the poten- 
tial, Vr(r), for the temperatures listed in Table U for the 
vector and pseudoscalar channels separately. In Eq. ([S]), 
standard symmetric lattice finite differences are used for 
the spatial and temporal derivatives. Figure [1] shows the 
central potential obtained for T/T c = 0.42 and 1.68 as a 
function of the time, r, appearing in Eq. (|5|). For each t 
value in FigUJ we have vertically shifted the data points 
so that V c (r/a s = 1) = 0. 

As can be seen, there is a good plateau where the po- 
tential Vc{r) is stable. The lack of a plateau at small 
times is presumably due to lattice artefacts caused by 
contact terms at the source. The upward trend of data 
points at large times and high temperature corresponds 
to time values close to the centre of the lattice which are 
contaiminatcd by backward moving states. We have con- 
firmed this interpretation by successfully modelling the 
effects of these backward moving states. 

The central values for the potentials are obtained from 
t = 6, 7, 7, 7 and 24 for N T = 20, 24, 28, 32 and 80 re- 
spectively. The resulting Vc and Vs are shown in Figs. 
[5] and [31 The left-hand error bars are statistical, and 
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the systematic uncertainty of choosing different values 
of t to define the potentials are depicted in the right- 
hand error bar. In Figj2]we include the Cornell potential, 
Vjr) = -f + £ + V with k = 0.52 and a = 2.34 GeV" 1 
[24| , as a point of reference. 
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FIG. 2: Spin independent (i.e. central) potential, Vc(r), for 
the temperatures in Table [T] obtained from eq.©. The data 
points have been shifted horizontally for clarity. The solid 
curve is a fit to the Cornell potential [24| (see text). 
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FIG. 3: Spin dependent potential, Vs(r), for the temperatures 
in Table U obtained from eq.©. The data points have been 
shifted horizontally for clarity. 

From Figl5]we see a clear temperature dependence. In 
the confined phase, T = 0.42T C , we see evidence of a 
linearly rising potential in agreement with the Cornell 
potential. As the temperature increases beyond T c , the 
potential flattens for large distances, in agreement with 
expectations of a deconfined phase. The spin dependent 
potential is plotted in FigJ3]and shows a repulsive core. 

We now compare our results with those using static 
quarks. There are two general approaches to extract 
the interquark potential between static (infinitely heavy) 



quarks, both of which have limitations. The first calcu- 
lates the free energy of a static quark pair as a function 
of their separation via various correlators of Polyakov 
loops [F3l - [l5| . However, the potentials thus derived suffer 
from either gauge dependence (in the case of the "singlet" 
channel), or do not reduce to the correct Debye screened 
potential in the perturbative limit (in the case of the "av- 
eraged" channel) [HI, H|| . The second uses Wilson loops 
or correlators of Wilson lines fl5l - fl7| and requires there 
to be good ground state dominance. However, this cre- 
ates tension because the temporal extent of the lattice, 
N T ~ 1/T is necessarily small at high temperature. As a 
result, precision results are difficult to obtain. 

The method discussed here is gauge invariant by con- 
struction and produces results with reasonable system- 
atics. Furthermore it calculates the potential between 
quarks with masses tuned to the physical charm. 

In Fig[J]we compare our results from Fig|3]with those 
obtained from static quark calculations - the singlet free 
energy [ijj] and the Wilson loop and line [T3|- We note 
a clear discrepancy between the our results and those 
obtained from static quarks. 
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FIG. 4: Comparison of Vc from this work with the singlet 
free energy calculation, Fi, from [l4| and the Wilson loop, 
Va, and Wilson line correlator, V\\, from [Vq . The error bars 
of the free energy data are smaller than the symbols. The 
data from this work follows the legend in Fig. [2] 

Conclusions - There is a significant body of theoret- 
ical work studying the interquark potential at non-zero 
temperature using model, perturbative and lattice (non- 
perturbative) approaches. The work outlined here uses 
a lattice simulation of QCD with two light dynamical 
flavours on an anisotropic lattice. We determine the char- 
monium potential at a variety of temperatures using rel- 
ativistic quarks tuned to the physical charm quark mass. 
This improves upon earlier lattice simulations performed 
in the static limit. It thus represents the first ab initio 
calculation of the charmonium potential of QCD at finite 
temperature. 

The method we use is based on the HAL QCD "time- 
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dependent" approach which obtains the real part of the 
potential from correlators of point-split operators fl9| . 
This allows the extraction of the potential without the 
need to first define the NBS wavefunctions by fitting the 
large time behaviour of the correlation functions. 

Our determination of the potential shows a linearly 
rising potential for T < T c and a clear temperature 
dependent flattening of the potential for T > T c . We 
demonstrate a significant deviation between our results 
and those obtained using static quarks via either the free 
energy or Wilson loops/lines. 

This work adds to previous charmonium studies per- 
formed by our collaboration with the same lattice param- 
eters [iH H(| and our earlier work on the potential using 
the HAL QCD wavefunction method [2l| . 

In forthcoming work we will simulate on significantly 
larger lattices with 2+1 light quark flavours. We also 
hope to extend our work to the potential between heavier 



quarks using the NRQCD approach [27i - f29j . 
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